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1. Introduction 

Higher-harmonic generation is a fundamental nonlinear 
optical process that has attracted intense interest ever 
since the beginnings of nonlinear optics [1, 2]. Genera¬ 
tion of optical waves via nonlinear interaction provides 
an effective approach to explore the properties of light 
matter interaction at fundamental level and, equally im¬ 
portant, has found key applications in many areas of 
science and technology. While most of the early stud¬ 
ies of higher-harmonic generation focused on wave in¬ 
teraction in homogeneous nonlinear optical media, re¬ 
cent advances in nanofabrication techniques and the ad¬ 
vent of metamaterials have considerably broadened the 
experimental and theoretical framework in which these 
nonlinear optical processes are explored. In particular, 
second-harmonic generation (SHG) from isotropic [3-5] 
and chiral [6, 7] metasurfaces, layered media coupled 
to arrays of plasmonic particles [8], and quantum engi¬ 
neered plasmonic metasurfaces [9] has been observed. 

A key prerequisite condition for achieving efficient 
SHG is that the wavevectors of the interacting waves 
are phase-matched. One particularly efficient method to 
phase match optical waves is by using diffraction grat¬ 
ings. In this approach, the residual wavevector imbal¬ 
ance is canceled by the grating wavevector, k = 27r/A, 
where A is the grating period. Importantly, the wide 
range of geometrical and materials parameters charac- 
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terizing a diffraction grating allows one to employ such 
optical devices to achieve wavevector phase matching in 
a multitude of optical wave configurations. More impor¬ 
tantly in this context, strong optical field enhancement 
and, consequently, a significantly more efficient nonlin¬ 
ear optical interaction can be achieved by using metallic 
gratings. Because generally there is no explicit analyt¬ 
ical solution to the problem of diffraction by periodic 
structures, efficient and accurate numerical methods are 
essential for the design of diffraction gratings with pre¬ 
defined spectral characteristics. Solving this problem in 
the nonlinear case is an even more daunting task, chiefly 
because the larger number of interacting waves and their 
intricate interplay significantly increases the complexity 
of the problem. 

The periodic nature of diffraction gratings renders 
frequency (wavevector) domain methods based on the 
Floquet-Bloch theorem to be particularly suitable algo¬ 
rithms for describing such optical structures. A common 
characteristic of these methods is the Fourier expansion 
of the optical field in a suitably chosen basis of modes, 
the main unknowns to be found being the corresponding 
Fourier coefficients. Among the most used methods in 
this class one should mention the Fourier modal method 
(FMM), also known as rigorous coupled-wave analysis 
(RGWA) [10-13, 15], the C-method [16, 17] for corru¬ 
gated gratings, and Green’s function type of methods 
[18-20]. On the other hand, there is a scarcity of numer¬ 
ical methods for SHG and other nonlinear interactions 
in periodic structures and the convergence and runtime 
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behavior of the few that do exist are not as well char¬ 
acterized. These methods are either problem specific 
algorithms [21] or extensions of linear methods to the 
nonlinear case, e.g. the FMM [22-24], finite-difference 
time-domain method [25, 26], Green’s function method 
[27], and multiple scattering matrix method [28, 29]. 

In this work we show how the linear generalized 
source method (GSM) [20, 31] can be adapted to de¬ 
scribe the SHG in diffraction gratings containing non- 
centrosymmetric materials. We choose the linear GSM 
as a starting platform for our new method because it has 
proven to be an efficient alternative to the FMM for the 
analysis of one- and two-dimensional (ID, 2D), corru¬ 
gated, periodic optical structures. Its optimal runtime 
is of order No\og{No), where No is the total number 
of diffraction orders, as compared to for the con¬ 
ventional FMM. The key ingredients needed to achieve 
this remarkable computational efficiency is the use of 
a matrix-multiplication based iterative solver, namely 
the generalized minimum residual (GMRES) method, 
to solve an underlying linear system of coupled inte¬ 
gral equations and a fast Toeplitz-matrix-vector multi¬ 
plication performed by means of a fast-Fourier transform 
(FFT). In addition, its implicit. Green’s function type 
formulation of Maxwell equations allows for a natural in¬ 
corporation of nonlinear source terms, which is the main 
challenge in computing the SHG response. 

The rest of the paper is organized as follows: in Sec¬ 
tion 2 we introduce the formulation of the nonlinear 
GSM, including the correct Fourier-factorization rule, 
and present the proper discretization of the underlying 
equations that allows one to use a fast iterative solver. 
Then, in Section 3, the convergence and runtime behav¬ 
ior of the method are studied. The method is applied to 
optimize the radiated second-harmonic (SH) by a diffrac¬ 
tion grating coupled to a slab waveguide in Section 4 
before final conclusions are drawn in Section 5. 

2. Derivation of the nonlinear generalized source 
method 

We start the derivation of the nonlinear GSM from the 
governing equations for SHG in the frequency domain 
and describe the general solution strategy. We then de¬ 
rive the nonlinear GSM for 2D periodic structures fol¬ 
lowing [20] and extend its scope to the SHG in optical 
gratings containing non-controsymmetric materials. 

2.A. Physical model for second-harmonic generation 

We consider a quadratically nonlinear optical medium 
and a time-harmonic electric field, £^(x, t), propagating 
in the medium at the fundamental frequency (FF), cJo, 

£:(x,t) =E(x)e^^°^-hc.c., (1) 

where x and t are the position vector and time, respec¬ 
tively. Through nonlinear interaction with the optical 
medium, this field gives rise to a nonlinear source polar¬ 
ization, P(x,t) = P^^(x) exp(if2t) + C.C., which oscil¬ 
lates at the SH frequency, Q = 2uJo. This polarization is 


the source of the electromagnetic field generated at the 
SH and is related to the excitation field via a third-order 
nonlinear susceptibility tensor, 

py(x) = (2) 

a, 13 

In this equation and what follows, Greek indices take 
the values of the Gartesian coordinates x, y, and z. 

Under the undepleted pump approximation, the ini¬ 
tially nonlinear problem in the time-domain is solved 
in three steps: i) Galculate the field at the FF, cjo; ii) 
evaluate the nonlinear polarization generated at the SH 
via Eq. (2); and in) calculate the field at the SH. The 
first and last steps are performed using the linear and 
nonlinear (extended) versions of the GSM, respectively. 

2.B. The GSM for inhomogeneous problems 

At both the EE and SH, the electromagnetic field is gov¬ 
erned by the electromagnetic wave equation with a phys¬ 
ical source current, J®^^(x), 

V X V X E — (3) 

where for convenience we have suppressed the spatial 
dependence of the variables and we introduced a generic 
frequency, cj, that takes the values cj = cjo (cj = U) at 
the EE (SH). In Eq. (3), ^(x) is the spatial distribution 
of the electric permittivity that defines the grating. 

The GSM is a rigorous method for solving Eq. (3). It is 
based on the decomposition of s{x.) defining the grating 
into a simple background structure with permittivity, 
56(x), and the difference structure, A5(x) = 5(x)—^^(x). 
The background has to be chosen such that the lin¬ 
ear solution operator characterizing the background 
structure problem, namely the operator that associates 
a source J with the corresponding solution of Eq. (3) 
with 5 = 56, is known. This means that 

E = N6(J) (4) 

solves Eq. (3) for = J and 5 = 56- 

In order to find the fields for an arbitrary permittivity 
5(x) one rewrites Eq. (3) as 

V X V X E - wVoebE = itofio + J®®"(E)] , (5) 

where the term J®®"(E) = —iuiAs'E, which depends on 
the solution E itself, is called generalized source. Using 
the solution operator this can be rewritten as 

E = N6 (J""^ + JS"^(E)) = N6 (J^°^(E)) , (6) 

with J^°^(E) = + Js®^(E). At the EE we choose the 

external source term such that the external field, 

Eext^ 

is an incident plane wave, although other choices 
are possible. At the SH the external source term is given 
in terms of the nonlinear polarization: 

So far, no particular assumptions about the structure 
under investigation have been made and thus the for¬ 
mulation (6) is valid for any structure. We consider a 
2D corrugated grating, as depicted in Eig. 1. It consists 
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Fig. 1. Setting for the GSM and closeup of the unit cell of 
a 2D corrugated grating containing nonlinear material. The 
grating is described by ^(x) and and is under plane 

wave incidence with wavevector ko. 


of the grating region 0 < ^ < /i, with grating vectors 
Ki and K 2 lying in the transverse {x^y) plane, and a 
given permittivity distribution, ^(x), which is a peri¬ 
odic function of the x and y coordinates. This slab is 
sandwiched in-between the cover and substrate, which 
consist of linear optical materials with permittivity 5c 
and 5s, respectively. We assume that only the grating 
region contains nonlinear material, i.e. = 0 if z < 0 
OT z > h. Due to its periodicity, the permittivity of the 
structure can be decomposed in a Fourier series. 


j(x) = y] £n{z)e 


i[{niKtx+n2K2x)x+{niKxy+n2K2y)y] 


where the sum over n = ( 711 , 722 ) is to be understood 
as a double infinite sum over the integers ni and 712 , 
711,2 ^ Because 5(x) is a periodic function in the 
transverse plane, one can use Bloch theorem to express 
the electric field as a Fourier series of functions that are 
pseudo periodic on the transverse coordinates x and y, 

00 

E(x)= (7) 

n= —00 


where k^^/y = k^^/y + niKi^/y + n 2 K 2 xiy are the pro¬ 
jections of the wavevector of the 71 ^^ order diffraction 
mode. The principal direction of the central diffraction 
order, which is defined by 711 =712 =0, is given by ko = 
/Cc(sin^cos(/), sin^sin^, cos^), where kc = cJoy/^cMo is 
the wavenumber in the cover region. For oblique inci¬ 
dence, this leads to a phase-shift of the electric field E(x) 
Qf g«(/coccAi+/coyA 2 ) where Ai and A 2 are 

the periods along the x and y axes, respectively. 

Since this analysis does not include nonlinear optical 
effects in the cover and substrate materials and since 
the nonlinear polarization P^^(x) arises due to a peri¬ 
odically distributed field (7) at the FF, both the external 


current at the SH and the generalized 

source term can be expressed as Fourier series. 


- iVt{e - Sb)E 

00 

= E (8) 


The coefficients depend on all field terms E^(^) 

- £b{z)] E(2:)}„ 

00 

= -iLO E - ^b{z)]^_^Em{z). 


m= — oo 


The solution operator N 5 for a periodic source cur¬ 
rent, in homogeneous space with con¬ 

stant permittivity 55 reads [18, 20 ]: 

En{z) = ^ J Gn{z,z')jn{z')dz', (9) 

with the tensor Green’s function 


G„(^, z') = (e+ef + h>E) Hiz-z')e<^^^-y^ 

-sf;^e^e^S{z - z'). (10) 

Here, and are TE and TM component vectors, 
respectively, of the ti^^ Fourier term of the electric field 
En, H{') denotes the Heaviside step function, and 
means the transpose of a vector a. The unit vector in z- 
direction is denoted by e^. Also, the dispersion relation 

l^nz — ~ {^nxY ~ {^nyf‘\ holds. OuC CaU SCO 

that the tensor gives rise to upwards and downwards 
propagating plane waves and a stationary term. In order 
to express the solution in terms of a superposition of 
plane waves, we eliminate the contribution of the (5-term 
in Eq. (10) by defining the modified electric field 

E^{z)=E^{z)-I^X°\z). ( 11 ) 

lUJSh 

The Eourier components of this modified field, E^(^) = 
Q^Sin{z)^ can be expressed via their TE/TM amplitudes, 
a.n{z) = {aYz),a-^{z),al^{z),al^{z)), and the matrix 
Qn = (^n 5 5 )• Additionally, Eq. (11) leads to 

a useful property of its 2 )-component: 

E,{z) = (12) 

The total source can be rewritten as 


Jtot _ jext ^ _ jext _ 

= -iLo(P^^ + eE - SbE) = -ioj (D - D;,) ,(13) 

where = e^E and = 0 if w = wq. So far, an 
infinite number of terms in the Eourier series was con¬ 
sidered. However, it is well known ([11, 12, 32]) that 
the factorization of products of periodic functions with 
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a finite number of terms has to be performed with care, 
namely the correct approach is dictated by the continu¬ 
ity properties of the factors in the product. To be more 
specific, products of functions with non-simultaneous 
discontinuities (e.g., the tangential component of the 
electric displacement field = eE\\ + can be fac¬ 
torized using Laurent’s product rule, whereas continuous 
products of functions with simultaneous discontinuities 
(the normal component D± = sE_\_ + is continuous 
at interfaces, whereas 5 and E± can be discontinuous) 
can be factorized using the inverse rule. In the context of 
numerical methods for periodic structures in the Fourier 
space, such as the FMM and GSM, the violation of the 
correct factorization rules causes slow convergence with 
respect to the number of retained diffraction orders. In 
contrast to previous work [33], where the correct factor¬ 
ization was only used in the linear part of the calcula¬ 
tions, we now show how the correct factorization should 
be implemented in the nonlinear part of the calculations. 

We now present the derivation of this modified rule 
and the impact of its correct application on the accuracy 
of the calculations is investigated in Section 3. Thus, 
to illustrate why a modified factorization rule for the 
inhomogeneous wave equation is necessary, we consider 
the example given in [32], adapted to our case: 



X 




kl < 7r/2, 
kl > 7i'/2, 


p{x) = 0.1c[l + sin 



c, 

0 , 


kl < n/2, 
\x\ > 7r/2, 


/ N pM 
d{x) = f{x)g{x) + p{x) 


r 6(1 - |a;|/7r), 

I 26(1 - \x\/tt), 


kl < n/2, 
\x\ > 7r/2, 


Here, a, 6 > 0 and c > 0 are arbitrary constants. One 
can easily see that the function d is continuous, whereas 
/ and g have concurrent discontinuities at \x\ =7:12. For 
c 7 ^ 0, p is also discontinuous at \x\ = 7r/2. Otherwise, 
Li’s example [32] is obtained. Figure 2(a) shows the 
functions for a = 6, 6 = 2, and c = 3. 

In order to obtain the vector, [d], of Fourier coeffi¬ 
cients of d, three different factorizations are compared, 
the results being depicted in Fig. 2(b). The product rule 
yields [d^] = |/][p] + [p], where |/] is the Toeplitz ma¬ 
trix of the Fourier coefficients of /. The reconstruction, 
R([d^]), of the function d shows strong oscillations at 
\x\ = 7r/2 around the expectedly continuous function d 
itself. The reason for this behavior is that the product 
rule propagates the overshoot from R([/]) and R([p]) 
to the reconstruction of the product R (|/][p])- In addi¬ 
tion, the oscillatory behavior of R ([p]) at |x| = 7r/2 adds 
to the variations of the reconstructed product, R([d^]). 

A direct application of the inverse rule for the prod¬ 
uct yields [d*] = |l//]~^[p] + [p]- The oscillations of 
the reconstruction R ([d*]) are now less pronounced but 
still consist of additive contributions from both terms 
R- (p-ffPlg]) and R([p]). 


Fig. 2. a) Discontinuous functions / (thick brown), g (dashed 
green), and p (broken blue lines) and continuous function d 
(solid black), for a = 6,5 = 2,c=3. b) Reconstruction 
of Fourier series with 200 components around discontinuity 
obtained by the product rule (broken brown lines), inverse 
rule (dashed green), and modified inverse rule (solid blue). 


We propose a modification of the product rule: we 
rewrite the continuous function d{x) = f{x)g{x)-\-p{x) = 
f{x)[g{x) Ep{x)/f{x)] as a product of discontinuous 
functions in the spaee-domain. The two factors have 
complementary discontinuities, whence the inverse rule 
applied to this product yields the correct factorization, 
[d"^] = {[g] + [p/f])' This is illustrated by the 

reconstruction R([d’^]), which shows no spurious oscil¬ 
lations at the points of discontinuity. 

Armed with this correct factorization rule, let us now 
consider Eq. (13). Thus, we find that in the context of 
the nonlinear GSM, d, /, p, and p correspond to the 
continuous normal component of the displacement field, 
the permittivity, 5, E±^ and respectively. The 
correct factorization of tangential and parallel compo¬ 
nents of D therefore reads: 


oNL 




[Dx] = [l/eF^ + [py/e]). 


(14a) 

(14b) 


We specify again that [V] denotes the vector of Nq = 
(2A^i + 1)(2A'2 + 1) Fourier coefficients of V. The trun¬ 
cation of Fourier terms with ni = — A^i,...,A^i and 
722 = • • • 5 5 with Ni and N 2 integers, corre¬ 

sponds to a rectangular truncation pattern in the Fourier 
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space. Problem specific, non-rectangular patterns could 
also be used and may lead to faster convergence in cer¬ 
tain cases (see [11] for a study of different truncation 
patterns in the FMM). Accordingly, the Fourier coeffi¬ 
cients of the total source current, defined 

by Eq. (8) are given by 



where I denotes the identity matrix of respective size. 

It remains to reformulate these equations, which are 
given in terms of the (A, ||)-components of E, in terms 
of the Cartesian components of E, the latter ones being 
the actual variables used in the GSM formalism. To this 
end, a local normal/tangential vector field is defined so 
that for any vectorial quantity v we can write: 


into the definition (11) of E to obtain 


[4] = -D iT,,[E,] + T,y[Ey]) + CE, + K], 


where C = A — DF^^. This yields [E] resolved in terms 
of [E], which in matrix form is expressed as 


f{E.] 

[Ey] 


I 0 0 

0 I 0 

c-'Dr,, c-'Dr.y c-i 



[E.] \ 

[Ey] 

[E.] ) 


(16) 


Plugging this formula back into Eq. (15) provides the 
total source current, which consists of the generalized 
source dependent on the unknown field E and a physical 
source whose origin is the known polarization: 


Vx\ ( 

t-y = B , 

Vz J \v4, J 

where the orthogonal transformation matrix 

( cos (j) sin ^ cos (j) cos ^ — sin (j) \ 
sin (j) sin -0 sin (j) cos -0 cos (j) 

cos^p —sin^p 0 J 

is a concatenation of the unit vectors in normal, and 
both tangential directions h, -0, and 0, respectively. Us¬ 
ing this transformation we get (analogously to Appendix 
A in [20]): 

ba] = y] {ASa/3 - [E/s] - iujeb[p^], (15) 


where 



= B diag ( 




( \fJe] \ 

4 


V 


J 


The matrices A and D are defined as 


A = 

5 

— 1 and D = 

5 


'£bl 


.^b_ 


.^b. 


. e 1 


and Fq,^ are the Fourier images of the entries of the 
matrix 


r(^) 


( cos^ (p sin^ A sin cp cos cp sin^ A cos cp sin A cos pj 
sin (p cos (p sin^ A sin^ (p sin^ pj sin cp sin pj cos A 
cos 0 sin A cos A sin 0 sin A cos A cos^A 


In agreement with the modified factorization rule, the 
Fourier series coefficients [p^ / 5 ] are obtained by numeri¬ 
cal integration of the space-domain quotient of the func¬ 
tions p^(x) and e{x). 

In order to eliminate E in favor of E we insert 


M 

iujSb 


D {r,,[E,] + r,y[Ey]) + (Dr,, - A) [e,] - \pp] 


M = + ( 17 ) 

where the modified polarization reads 

K] = -[|^] + (A 5 „,-Dr«,)C-i[|^] ( 18 ) 

and the matrix W( 2 )) is defined as 

w = j -h OTxx ~ ^ DF^^C ^ j 

V -c-'Dr,, -c-iDr,y c-i-i J 

Here, Nq,^ = DFc^^C^^DF^^ combines the geometrical 
normal field information contained in F with the physical 
structure given by 5 . 

Combining the main results, namely Eq. (17) that ex¬ 
presses the source in terms of E and the Green’s function 
(9), yields a system of Eredholm integral equations of the 
second kind for the unknown amplitudes, 

nh 

a.n{z) = a°^{z) + Rn(^,^')EQ:a,n 

No 

X ^ ,nm (.z')Q/S:,m^m{z')dz', ( 19 ) 

m= — No P 

for n = 1,..., No, where the sum over m = (mi, m 2 ) is 
to be understood as a double finite sum over the inte¬ 
gers mi = -Ai,..., Ai and m 2 = -A 2 ,..., ^ 2 ^ 
and Q/3:^rn denote the a-column and P-iow of and 
Q^, respectively. The matrix Rn{z,z') G incor¬ 

porates the propagation of TE- and TM-polarized, up¬ 
ward and downward plane wave amplitudes of the 
Eourier component inside the slab-background structure 
and their reflections at the top and bottom interfaces 
[see Eqs. (61) and (62) in [20] for the exact formula]. 
The known values of 8i^{z) are either the amplitudes of 
the incident plane wave electric field in the background 
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structure at the FF or the amplitudes of the external 
source field at the SH; the latter are given by 

W = / Rn{z,z')QlnpPrii.^')dz'. (20) 

Jz'=0 

One can easily see that the modified factorization rule 
does not affect the generalized source but only changes 
the external polarization. This means, that all algorith¬ 
mic enhancements, that make the linear GSM an effi¬ 
cient numerical method, can be applied here as well. 

Using the midpoint rule, the system (19) of integral 
equations is discretized at Ni equidistant points, Zp = 
{p — 0.b)h/Ni^ p = 1,..., This results in a system of 
linear equations. 


Ni 


a.np — ^ ^q) Q:c 


q=l 


No 


^ ^ ^ ^ al3,nm{^q )^l3:,m^mqi ( 21 ) 

m=—No /3 


for n = 1,..., Nq^p = 1,..., Ni. The total number of 
unknowns in Eq. (21) is and usually becomes pro¬ 

hibitively large for direct solvers to be used in practical 
situations of interest. 

In order to take full advantage of an iterative solver, 
the system (21) is reformulated to (see also [20]): 

a= (I + RQUA-^Q) a^ (22) 


where a and a^ are the vectors containing the dis¬ 
crete unknowns, a^p = 3in{zp), and the known ampli¬ 
tudes, a^p = a^(zp), respectively, A = M — QRQU, 
and Q and Q are block-diagonal matrices with and 
on their block-diagonal, respectively. The matrix 
R G C^NoNixANoNi corresponds to Rn(^p,^g)- Its (4(n — 

l)7V,+4(p-l) + [l,..,4], 4(n-l)7V,+4(g-l) + [l,..,4])- 
entries are given by Rn(^p, Zq). Hence, R is diagonal with 
respect to the Fourier component n and has Toeplitz- 
structure with respect to the layer indices p and q. The 
matrices 


M = 


(iTlb] 

0 

0 ^ 


0 


raid 

0 

sin^ ^ + cos^ -0, 



0 1 

[i]Kh 


-[?1 

£ 

, G = l- 

e 

p6] 
1 s 1 

5 


U = 


f AMcca, -h G|^|ra, 

Fr., 


AM^/y + G 

FF, 


J- xy 


\ 


GF, 

GF, 

~ l^iJ 


are defined such that W = UM“^ but they do not 
contain inverted submatrices themselves. Additionally, 
U and M are block-Toeplitz-Toeplitz-matrices with re¬ 
spect to the index pair n = (ni, 77 , 2 ). Due to the Toeplitz- 
property of its submatrices, multiplication of a vector b 


of size SNqNi with the SNqNi x SNqNi system matrix A 
can be performed in 0{NoNi log(A^o) log(A'/)) [30] opera¬ 
tions instead of 0{N‘^Nf) for the standard matrix-vector 
multiplication. Therefore, the use of an iterative solver 
based on matrix-vector multiplications like GMRES is 
highly beneficial. Although the reformulation Eq. (22) 
removes all inverted submatrices from the system matrix 
A itself, the evaluation of the modified source polariza¬ 
tion [p^] in Eq. (18) still requires a one-time inversion 
of matrix C. In practice, the computational cost of this 
nonrecurring inversion is, however, negligible compared 
to the overall runtime of the algorithm, which is domi¬ 
nated by the inversion-free iterative solution process. 

3. Numerical analysis and convergence studies 

In this section we examine the convergence and runtime 
properties of the linear and nonlinear parts of the GSM 
using cylindrical gratings as test examples 2D made ei¬ 
ther of dielectric or metallic materials. Importantly, we 
illustrate the impact of the application of the correct 
factorization rule as compared to the old factorization. 

The unit cell of the grating under consideration is as¬ 
sumed to be rectangular and, for simplicity, the two pe¬ 
riods are set to be equal, Ai = A 2 = A = 2 pm. It 
contains a cylindrical disk of height, h = 200 nm, and ra¬ 
dius, r = 0.25A = 500 nm. The disk and substrate have 
refractive index, = 2, whereas the cover region is vac¬ 
uum with nc = 1. We assume an isotropic nonlinear sus¬ 
ceptibility tensor, xfki = X^jkhh with y = 10“^ m 
The incident plane wave at the fundamental wavelength 
of Aff = 1.5 pm, is specified by the angles 0 = 30° and 
<p = 70°, and a 7 /^ = 45° polarization. 

The optical response of the device at the EE and SH 
is computed by the linear and nonlinear GSM, respec¬ 
tively. We perform simulations for an increasing num¬ 
ber of diffraction orders No = (2A^i + 1 ){ 2 N 2 + 1) for 
Ni = N 2 ^ {3,5, ...,23} and increasing number of lay¬ 
ers, Ni G {12, 25, 50,..., 400}. In order to assess the 
efficiency of the algorithm, we compute the following 
quantities: the reflection and transmission coefficients, 
rFF/sh rj^FF/SH^ respectively, at the EE/SH, the 

number and of GMRES iterations necessary 
to solve Eq. (22) at the EE and SH, respectively, and 
the GSM runtime needed to complete the linear (t^^) 
and nonlinear {t^^) part of the computations. 

Eor the validation of the linear calculations results, 
they are compared to data obtained by using an effi¬ 
cient normal vector field formulation [13] of the EMM, 
which is also commercially available [14]. Eor the pre¬ 
sented comparison, both methods were implemented in 
MATLAB and all simulations were performed on a single 
2.5 GHz processor to allow a fair comparison of simula¬ 
tion times. 

The results of this analysis are summarized in Eig. 3. 
The dependence of the transmission coefficient at the 
EE, on the number of diffraction orders, computed 

by using the GSM for several values of the numbers of 
layers Ni in which the grating is divided, is presented 
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Fig. 3. a) Transmission efficiency obtained using GSM for in¬ 
creasing Ni [colors, legend from b)] and FMM (black). Inset: 
discretization error of GSM vs. No. b) Runtime of GSM for 
increasing Ni (colors) and of FMM (black), c) Gonvergence 
behavior of radiated SH in the direction of the incoming plane 
wave using the two different factorization rules. 


in Fig. 3(a). For comparison, we present the same de¬ 
pendence determined by using the FMM. Both methods 
show a similar convergence pattern with respect to Nq 
and the computed transmission coefficient rapidly ap¬ 
proaches the asymptotic value of = 0.917772. For 
the investigated lossless device the GSM conserves the 
incoming optical power up to a relative error of 10“^ for 
low No < 1000 and the power conservation improves to 
10“^ as No further increases. 

The inset in Fig. 3(a) illustrates the error in the trans¬ 
mission due to the discretization of the integral equation, 
defined as AT = calculated for several 

values of orders, No G {49, ...,2209}, and for increasing 
Ni. Given sufficiently many orders, e.g. No = 2209, 
quadratic decrease of AT is achieved, which is to be ex¬ 
pected since using the midpoint rule in the discretization 
of the integral equation implies quadratic convergence. 
For smaller e.g. No = 225, the convergence plateaus 
at a small number of layers. This is because the re¬ 
formulated linear system (22), which enables the use of 
the accelerated matrix-vector multiplication, is only No- 


Table 1. Gharacteristics of the nonlinear GSM for dielectric 
and lossy materials 


n 


A?^ 

^ iter 

rjnF F 


rj-\SH 

1.45 

21 

38 

0.972998 

2.6 • 10“® 

2.88 • 10“^’’ 

2 

39 

100 

0.917772 

8.0 • 10-'^ 

1.16- 10“^’’ 

3 

105 

400 

0.825020 

5.8 • 10“® 

4.57 • 10“^® 

1.5 + lz 

96 

157 

0.620149 

0.283418 

7.41 • 10’^® 

0.36+6.5Z 

2160 

1968 

0.381113 

0.193255 

5.99 • 10“^® 


asymptotically equivalent to the original system (21). 

Figure 3(b) depicts the runtime necessary to com¬ 
pute the solution of the problem at the FF for increas¬ 
ing No and for different A/. The linear GSM exhibits 
nearly optimal runtime as increases log-linearly with 
both No and A/. This is because is chiefly deter¬ 
mined by the computational work needed to iteratively 
find the solution of Eq. (22) using GMRES, which in 
turn depends on the number of iterations A^^^ and, if 
« NM, on the time required for the system- 
matrix vector multiplication; this latter quantity is of 
order 0{NoNi log(Ao) log(A/)). Eor the investigated de¬ 
vice, only = 39 GMRES iterations are necessary 

to solve Eq. (22), where is nearly independent of 

No and Ni (see also Table 1, second row). The black line 
in Eig. 3(b) indicates the runtime of the EMM, 
which follows a 0{N^) dependence. This is the expected 
behavior as the EMM involves finding the solution of an 
eigenproblem of size A^. The asymptotic advantage of 
the GSM becomes evident, e.g. at Ni = 200, where for 
No = 1225, becomes larger than . 

Eor the computations at the SH, the GSM compu¬ 
tational time still follows the ©(A^A/log(Ao) log(A/)) 
trend; however, as compared to the EE, approximately 
3x more iterations, A^f^ = 100, are necessary to solve 
Eq. (22) at the SH. We identify two possible reasons for 
this behavior: Eirst, we found that for plane wave ex¬ 
citation, the smaller the wavelength the more iterations 
were necessary for GMRES to solve Eq. (22) and second, 
both the r.h.s. of this equation and its solution are more 
inhomogeneous at the SH as compared to the EE. 

The results for the generated SH, computed with 
Ni = 400, are depicted in Eig. 3(c). Specifically, this 
figure shows the radiated power at the SH in the direc¬ 
tion of the incoming wave, normalized to the incident 
power at the EE. The data corresponding to the solid 
line was obtained by using the modified inverse rule for 
the source polarization and shows remarkably fast con¬ 
vergence to a value of = 1.163-10“^^ with increasing 
No. By contrast, when using the unmodified rule 
does not reach its asymptotic value even for as many as 
No = 2209 orders, despite the fact that the EE solution 
is already converged if No ^ 750 [see Eig. 3(a)]. Note 
that the old and correct factorizations tend to slightly 
different values. Such a behavior has already been re¬ 
ported for linear calculations in the EMM ([11]). 

Einally, we investigated the convergence properties of 

































the GSM when it was applied to simulate diffraction 
gratings made of three different lossless dielectric mate¬ 
rials, a lossy dielectric, and metallic components. The 
results, shown in Table 1, provide valuable insights into 
the characteristics of the GSM. Thus, clear trends can 
be observed in the case of lossless dielectric gratings: the 
higher the refractive index n is the more GMRES itera¬ 
tions are necessary at the FF and the number of it¬ 
erations at the SH increases faster than . The 
numerical absorption, = 1 — ^ which in 

the case of gratings made of lossless materials quantifies 
the degree to which the algorithm conserves the energy, 
decreases with the index of refraction. The lossy mate¬ 
rial with n = 1.5 + i can be efficiently analyzed by the 
GSM: the is relatively small and increases only by 
50% for the calculation of the SH electromagnetic field. 

The metallic grating requires a more detailed discus¬ 
sion. The derivation of the GSM formally allows com¬ 
plex values of the permittivity e{x.) and the algorithm’s 
implementation does show that the method can be ap¬ 
plied to simulate diffraction by lossy structures (see the 
n = 1.5 + i example in Table 1). However, we found out 
that the very large number of required GMRES itera¬ 
tions makes the method particularly slow when applied 
to structures containing very lossy components, namely 
a material with n = 0.36 -1- 6.5i. We also found out 
that for metals A^iter increases strongly when Nq in¬ 
creases (the presented results for n = 0.36 + 6.5i are 
for No = 529; increasing Nq would require A^iter > 5000 
to achieve convergence of GMRES). Summing up these 
results, the GSM in the presented form is most efficient 
and most suitable for shallow, low refractive index con¬ 
trast, dielectric devices. However, an improvement of 
the GSM using curvilinear coordinates [34] allows one 
to efficiently calculate metallic structures as well but we 
have not implemented this version of the algorithm in 
the nonlinear case. 

4. Application to SHG in textured slab waveguides 

In this section we illustrate the versatility of our method 
by showing how it can be applied to a problem of practi¬ 
cal importance, namely to the maximization of the SHG 
in a phase-matched configuration. Because of the weak 
nature of nonlinear optical interactions, the radiated SH 
in practical settings is orders of magnitude smaller than 
the linear excitation. Therefore field enhancing mecha¬ 
nisms are particularly important in this context as they 
can be employed to increase the efficiency of SHG in 
practical applications. To this end, we consider the 
SHG in a diffraction grating resonantly coupled to a slab 
waveguide, both optical devices being made of quadrat- 
ically nonlinear optical materials. 

To add specificity to our problem, we assume that the 
nonlinear material in the device is GaAs (refractive in¬ 
dex n = 3.4), which due to its large optical nonlinearity 
is a particularly suitable choice for our device applica¬ 
tion. It crystallizes in the zincblende lattice type, so 
that it belongs to the crystal point group 43m. Hence, 
the non-vanishing components of its second-order sus¬ 


ceptibility tensor are Xxyz = Xxyz = Xyxz — Xyzx — 
X^zxy = X^zyx = 7.4 X [35], which is a rel¬ 

atively large value as compared to that of most other 
non-centrosymmetric materials. As substrate we choose 
GaF 2 {us = 1.4) and consider air {uc = 1) as cover, 
which provides the high refractive index contrast needed 
to achieve good field confinement in the slab waveguide. 

In order for the incident plane waves to couple and 
interact with the slab waveguide modes a periodic op¬ 
tical grating that breaks the translational symmetry of 
the structure is placed on top of the waveguide. As a 
result the propagation constant, /3jy(A), of the z/th slab 
waveguide mode is folded within the first Brillouin zone 
of the grating. When the difference between the tan¬ 
gential component of the incident wavevector, ky, and 
the propagation constant /3jy(A) is equal to a multiple 
of the grating vector this waveguide mode can be res¬ 
onantly excited. In particular, one expects to observe 
strong field enhancement in the slab waveguide at the 
corresponding resonance wavelengths and consequently 
significantly stronger nonlinear optical interaction, an 
effect that has been recently used to demonstrate en¬ 
hanced optical absorption in such optical devices [36]. 

In the context of SHG there are two kinds of reso¬ 
nances that can affect the nonlinear optical response of 
the device: linear resonances correspond to the grating- 
induced, resonant excitation of a waveguide mode at the 
FF by a plane wave incident at a certain angle and wave¬ 
length, Xff- We expect then that the optical field at the 
FF is strongly enhanced and confined inside the waveg¬ 
uide, thus a stronger nonlinear polarization is created, 
which in turn excites a stronger field at the SH wave¬ 
length, XsH = Xff/‘^- We call this phenomenon of di¬ 
rect enhancement of the nonlinear response at the SH 
due to the excitation of a linear resonance at FF an in¬ 
herited resonance at Xsh- The second type of nonlinear 
resonances are intrinsic resonances^ which are observed 
when a waveguide mode is excited at the SH wavelength 
by the dipoles associated to the nonlinear polarization. 
In this case the field at the SH is confined inside the 
waveguide and therefore shows remarkably large radi¬ 
ated SH intensity. Note that the intrinsic resonances 
are not necessarily grating coupled with the radiative 
continuum and therefore they can be viewed as nonlin¬ 
ear dark modes of the system [37]. In particular, one 



Fig. 4. a) Unit cell of a one-dimensional binary grating and 
b) two-dimensional cylindrical grating on top of a slab waveg¬ 
uide made of nonlinear GaAs on a CaF 2 substrate. 
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expects that these resonances have large quality factor. 
We call the inherited and intrinsic resonances nonlin¬ 
ear resonances, since they are observed at the SH. For 
the sake of simplicity, the modes causing the linear and 
nonlinear resonances are denoted as linear and nonlin¬ 
ear modes. The SH wavelengths Xsh for which both 
types of nonlinear modes are excited simultaneously are 
of particular interest because at these wavelengths, this 
doubly resonant mechanism of SHG, can lead to a re¬ 
markably large nonlinear optical response of the device. 
In the remaining of this section, these ideas are illus¬ 
trated by studying the nonlinear optical response of ID 
and 2D implementations of these devices. 

4.A. One-dimensional binary gratings 

We first investigate a ID grating, which although having 
a simple structure, provides us the necessary insights 
to completely understand the more complex 2D design. 
Thus, let us consider a ID binary GaAs grating on top 
of a GaAs slab waveguide, placed onto a GaF 2 substrate, 
see Fig. 4(a). The height of the slab waveguide, is a 
free design parameter and varies from 0.25 pm to 0.5 pm. 
The period of the grating with filling factor p = 0.25 is 
A = 2 pm, the height of the grating region is chosen in 
relation to the height of the slab waveguide, hg = 

At the FF, we consider a TM polarized plane wave 
with wavelength XpF ranging from 3.6 pm to 4.4 pm un¬ 
der 0 = 30° incidence. By considering the slab waveg¬ 
uide region as homogeneous region of the periodic grat¬ 
ing, its nonlinearity is also calculated by the nonlin¬ 
ear GSM. For the calculation of the presented results 
No = 21 diffraction orders and Ni = 138,..., 275 GSM- 
layers were used, each with height of about 2 nm. 

Figure 5(a) shows the transmitted intensity at the FF 
for a TM polarized incident plane wave with wavelength 
XpF’ The spectrum shows a smooth dependence ohXff 
and except for a narrow stripe that corresponds to 
a sharp decrease of the transmission. This change in 
transmission can clearly be attributed to the excitation 
of a TMo waveguide mode. Thus, an analytically cal¬ 
culated resonance wavelength [33], which is represented 
by the pink line in the spectrum in Fig. 5(a), is in close 
proximity to the band of reduced transmission seen in 
the spectrum. The analytical calculation is obtained for 
a infinitesimally small perturbation of the slab waveg¬ 
uide, whence the agreement is nearly perfect for low 
hw and slowly deteriorates as and therefore hg in¬ 
creases. Moreover, Fig. 6(a) shows the ^-component of 
the magnetic field at the FF, inside and around 

the textured waveguide, which is indicated by blue lines, 
corresponding to the black cross (1) in Fig. 5(a), i.e. to 
Xff = 3.861pm and = 345 nm. This field profile 
clearly matches that of the TMq mode. Note that for a 
ID grating under non-conical incidence, the TE and TM 
components of the electromagnetic field decouple. Ac¬ 
cordingly, there is only a single TM mode excited under 
TM-incidence, in the considered H^-Xff domain. 

In order to better visualize the SHG enhancement. 



AFp(/im) 



1.80 1.93 2.07 2.20 

AsH(/um) 



Fig. 5. Spectra of transmission coefficient for the binary tex¬ 
ture determined for different waveguide heights under nor¬ 
mal incidence and for TM polarization, a) Spectra at the 
fundamental wavelength (pink dispersion curve indicates an¬ 
alytical calculation of TMq mode), b) Spectra at second 
harmonic wavelengths (logarithmic scale), c) Zoom-in of the 
region around a point where an intrinsic and inherited SH 
resonance dispersion curve cross. 

we plot in Fig. 5(b) the logarithm = logiQT^^ of 
the SH intensity, emitted in the direction of the inci¬ 
dent wave, for the same incident fundamental field as in 
Fig. 5(a). The transmission coefficient, is normal¬ 

ized to the SHG in a bulk layer of GaAs having the same 
volume as the textured grating. This figure clearly shows 
an increase of the SH intensity due to the excitation of 
the inherited TM modes and intrinsic TE modes. No 
intrinsic TM modes are excited in this case. This is due 
to the particular symmetry properties of for a TM 
polarized fundamental field E(x) = (F^a,(x), 0, E;2(x)), 
only the ^-component (i.e. the TE component) of the 
nonlinear polarization is non-vanishing. 

At Xsh = 2.196 pm and = 460 nm, corresponding 
to the blue cross (3) in Eig. 5(b), one can achieve an 
intensity enhancement of 3 due to the 

excitation of an intrinsic TEi mode. Although the in¬ 
tensity is only three times higher than in the case of the 
bulk material, it is still significantly larger than the sur¬ 
rounding designs of textured waveguides. The excitation 
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Fig. 6. Magnetic and electric field distributions at simultane¬ 
ous resonance wavelength Xff = 3.45 pim for hw = 381 nm. 
a) \Hy^\ of the fundamental TMo-resonance at the FF. b) 
\Ey^\ of simultaneous resonance at SH. c) \Ey^\ of intrinsic 
TEi mode at the SH wavelength Xsh = Aff/2 = 1.725 pm 
for a TE polarized incident plane wave. 


of the inherited mode yields much stronger enhancement 
of SHG, namely at Xsh = 2.011pm and = 451 nm 
(red cross (2)), = 10^-^^. 

The strongest SH enhancement of = 10^’^ is 

achieved for simultaneous excitation of the inherited 
TMo mode and intrinsic TEi mode at Xsh = 1.931pm 
and hw = 345 nm (black cross (1)). A zoom-in of the 
region of the doubly resonant excitation is presented in 
Fig. 5(c). The electric field profile Ey^ at the SH for 
this particular configuration, plotted in Fig. 6(b), shows 
high field enhancement inside the waveguide-grating re¬ 
gion with two pronounced maxima at the top and bot¬ 
tom of the structure. This field profile can be ex¬ 
plained by inspecting the distribution of the field Ey^ 
of the intrinsic TEi mode at Xsh = 1.725 pm and 
hw = 381 nm. Thus, for TE-polarized incident waves 
with Xff = ‘^^sh = 3.45 pm, the electric field around 
the textured waveguide closely resembles the profile of 
a TEi mode, as can be seen in Fig. 6(c). However, the 
field enhancement is not as strong as in the case of the 
linear TMq mode, as per Fig. 6(b). 

This analysis fully explains the dominant feature of 
the SH response of the ID device when the two types of 
resonances are simultaneously excited, namely the over¬ 
all SH intensity is mainly determined by the strong en¬ 


hancement inside the waveguide of the fundamental field 
due to the linear resonance. The field profile, however, 
is determined by the profile of the intrinsic TEi mode. 

4.B. Two-dimensional cylindrical gratings 

The ideas presented in the previous section can be easily 
extended to a similar 2D device. To illustrate this, we 
consider the GaAs grating in Fig. 4(b) with unit cell con¬ 
sisting of a cylinder with radius, r = 500 nm and relative 
height hg = O.lhw^ placed at the center of a square with 
side-length A = 2 pm. The GaAs waveguide is placed on 
a GaF 2 substrate and its height, hw^ varies from 660 nm 
to 802 nm. For the presented computationally expen¬ 
sive parameter sweep, a total of A^o = 81 diffraction or¬ 
ders and Ni = 350,..., 450 GSM-layers with individual 
height of 2 nm were used. By checking the convergence 
with respect to Nq for a fixed hw = 700 nm and vary¬ 
ing wavelength, we found that No = SI already yields 
qualitatively correct results. 

The results for normally incident, TM-polarized plane 
waves with wavelength Xff = 4.12 pm to 4.44 pm are 
summarized in Fig. 7. As in the ID-case, the transmis¬ 
sion spectrum at the FF reveals the resonant excitation 
of a TMo mode, seen in Fig. 7(a) as a dark stripe. 

The spectrum of radiated SH in the direction of the 
incoming wave, presented in Fig. 7(b), exhibits dis¬ 
tinct maxima along the dispersion curves of the in¬ 
herited and intrinsic resonances. For example, for 
the inherited TEo-resonance at Xsh = 2.126 pm and 
hw = 780 nm (blue cross (3)) the SHG enhancement as 
compared to a uniform GaAs layer with the same vol¬ 
ume is = 10^’^^ ~ 5.9. The inherited resonance at 



2.06 2.11 2.17 2.22 

AsH(/^m) 


Eig. 7. Spectra of transmission coefficient for the cylindri¬ 
cal texture determined for different waveguide heights under 
normal incidence, a) Spectra at the fundamental wavelength, 
b) Spectra at SH wavelengths (logarithmic scale). 
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-A/2 0 A/2 ^ 

Fig. 8. Magnetic and electric field distributions at simultane¬ 
ous resonance wavelength Xff = 4.213 pm and hw = 696 nm. 
a) \Hy^\ of the fundamental TMo-resonance at the FF. b) 
\Ey^\ of simultaneous resonance at SH. c) \Ey^\ of intrinsic 
TEo mode at SH wavelength Xsh = Aff/ 2 = 2.107pm for 
TM polarized incident plane wave. 


^SH = 2.201pm and = 763 nm (red cross (2)) yields 

rj^SH _ 2q3.95 

We have designed a 2D periodic device in which 
inherited and intrinsic resonances can be simulta¬ 
neously excited, the corresponding parameters being 
Xff = 4.213 pm and hyj = 696 nm (black cross (1)). 
This double resonance leads to an intensity enhancement 
of « 10^ as compared to the reference uniform layer. 
The spatial profile of the magnetic field, calculated in 
the x- 2 ;-plane passing through the center of the cylinder 
{y = 0) is displayed in Fig. 8(a). The field distribution 
closely resembles the profile of the TMo-mode, showing 
strong field enhancement and increased field confined in 
the waveguide as compared to the ID-case. Our simula¬ 
tions suggest that normal incidence and a larger waveg¬ 
uide thickness favor stronger confinement. Moreover, 
the intrinsic TEq mode profile is given in Fig. 8(c). As 
in the ID-case, it strongly influences the shape of the 
SH-field inside the textured waveguide, as per Fig. 8(b). 

5. Conclusion 

We have presented an extension of the generalized source 
method to analyze SHG in 2D periodic structures con¬ 
taining non-centrosymmetric, quadratically nonlinear 
optical materials. The linear and nonlinear parts of the 
calculations are decoupled in the undepleted pump ap¬ 


proximation and a three-step computational process was 
derived. We introduced the correct Fourier-factorization 
rule for inhomogeneous problems and its beneficial ef¬ 
fect on the convergence in the nonlinear calculations was 
clearly shown. In a benchmark structure the asymp¬ 
totic runtime estimate of 0(A^log(A^)A'/log(A^/)) was 
validated, which renders the GSM especially suitable for 
modeling low permittivity contrast, dielectric gratings. 

As a practical application we optimized the design of 
ID and 2D textured GaAs slab waveguides to achieve 
maximal radiated SH by simultaneous excitation of lin¬ 
ear and nonlinear slab waveguide modes. In a addition 
to an analytical estimate for the resonance frequencies 
of the textured slab waveguide, investigation of the field 
profiles allowed us to identify the nature of the reso¬ 
nances responsible for the remarkably high intensity en¬ 
hancement of more than 8 orders of magnitude in com¬ 
parison to the SH radiated by a uniform slab. Impor¬ 
tantly, the formalism presented in this study can be ex¬ 
tended to several other key nonlinear optical processes, 
including surface SHG from centrosymmetric materials 
and higher-harmonic generation. Our method can be 
further extended beyond the undepleted pump approx¬ 
imation by allowing a self-consistent coupling between 
the FF and SH optical waves. In a similar manner, Kerr- 
nonlinear optical materials can be considered as well, a 
topic that we plan to address in a future study. 
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